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Abstract 



We discuss the effect of chiral anomaly as a possible mechanism for triggering 
formation of domains of disoriented chiral condensate (DCC) in relativistic 
heavy ion collisions. The anomalous ir° — > 27 coupling and the strong, Lorentz 
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contracted electromagnetic fields of the heavy ions combine to produce the 
"anomaly kick" to the field configuration of the neutral pion field. We im- 
plement the effect of anomaly kick in our numerical simulation of the linear 
sigma model in a schematic way which preserves its characteristic features: 
the effect is coherent over a large region of space but is opposite in sign above 
and below the ion scattering plane. We demonstrate by detailed simulations 
with longitudinal expansion that the DCC domain formation is dramatically 
enhanced by the anomaly kick in spite of its small absolute magnitude. We 
examine the behavior of various physical quantities such as pion fields, the 
axial vector currents, and their correlation functions. Our results also provide 
useful insight into the mechanism and properties of DCC domain formation, 
in general. Finally, we discuss some experimental observables which can signal 

the anomaly induced formation of DCC. 
PACS: 25.75.-q,12.38Mh,11.30.Rd 
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I. INTRODUCTION 



One expects, on account of universality arguments, that quantum chromodynamics 
(QCD) with two massless quark flavors exhibits a second-order phase transition between 
a low-temperature phase, which shows spontaneous chiral symmetry breaking, and a chi- 
rally symmetric high-temperature phase [|T],0. The fate of this transition for the physical 
values of the u- and (i-quark masses and a third semi-light quark flavor (s), is currently 
under intense investigation by means of numerical simulations of the lattice gauge theory If 
the transition exists in the real world, e.g., in intermediate states in heavy ion collisions, an 
interesting phenomenon can occur: Large coherent domains of pion fields form due to the 
long range correlations associated with the second order phase transition. Unfortunately, it 
is unlikely that the second order transition persists with finite quark masses, as it is easily 
destroyed by a weak external magnetic field in the case of magnetization, and we do not 
expect that large domains form in an equilibrium situation 0. 

However, it was argued by Rajagopal and Wilczek |J that domains may be formed 
in energetic collisions where hot regions experience subsequent non-equilibrium evolution. 
They proposed an idealized quench approximation to model this non-equilibrium scenario [[3] . 
It provides a concrete realization for the formation mechanism of large chirally misaligned 
domains, the disoriented chiral condensates (DCCs). Such chirally misaligned coherent pion 
field domains have been discussed by many authors Q after the pioneering works f5)-|7| in 
the context of large neutral-charged pion fluctuations (Centauro events) that may have been 
seen in cosmic ray experiments |§. 

The central question in the field, i.e., whether the DCC forms in high energy hadronic 
collisions, has been discussed extensively. We summarize these discussions below, and then 
we point out that the chiral anomaly combined with the environment of relativistic heavy 
ion collisions greatly enhances the possibility of formation of DCC domains. 

The possibility of the formation of large DCC domains is investigated by means of nu- 



merical simulation PJ9|,|10|j using the linear sigma model [11] as a low energy effective theory 
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of QCD. Among the most elaborate of these simulations, one with longitudinal expansion 
|I(J indicates that large spatial domains of DCC can, indeed, develop for appropriate initial 
conditions. 

Disoriented chiral domains can only form when a region of space is cooled down rapidly 
and chiral symmetry gets spontaneously broken sufficiently fast, so that the chiral order 
parameter (qq) cannot adiabatically follow the shifting minimum of the effective potential. 
This is the reason for the quench scenario. To realize it, the hot debris formed in high 
energy hadronic collisions must somehow be cooled down quite rapidly. A one-dimensional 
expansion does not appear to be sufficiently fast [fL0yi2| , as is expected from the early 
calculation using the hydrodynamical model by Bjorken [|13j. Therefore, a rapid three- 



dimensional expansion would be required flT^H] to drive the chiral field far enough out of 
equilibrium. 

If the relics of such domains were observed in experiments where matter is heated above 
T c , e.g., in relativistic heavy ion collisions, they would provide evidence for the existence for 
the chiral phase transition in QCD in non-equilibrium environments. Because the domains 
decay into coherent multi-pion states, they are predicted to reveal themselves in highly 
characteristic pionic observables. The unusual distribution of the pion charge ratio R = 
N w o/N n has been discussed extensively in the literature |3|,^H7|JT5|| , and it is now the subject of 
active experimental investigations ||16| , |17| . However, it has been shown that the charge ratio 
is useful as a sensitive signal for DCC formation only if a very small number of independent 
domains are formed [fL8yi9| , indicating the need for more advanced analyses. The possibilities 
include the wavelet analysis |19[ and multiparticle correlations [p!5|,|16[ . If the DCC can be 
described by a squeezed state |2jJ of the pion fields, as expected by the parametric resonance 



mechanism pl]| , some peculiar two pion correlations are expected ||22|| . 

We demonstrate in this paper that the Adler-Bell-Jackiw chiral U(l) anomaly [23] ex- 



pressed as the Wess-Zumino term ]2j] in the linear sigma model does affect the chiral orien- 
tation of pion fields, and thereby enhances the DCC domain formation. As will be explained 
in the following sections, we formulate the effect of the anomaly as a "kick" to the neutral 



pion fields imparted by the electromagnetic fields of colliding relativistic heavy ions p5| , and 
then implement it into a numerical code simulating the linear sigma model |]IU|| . We show 
by examining the various physical quantities such as pion fields, the axial vector currents, 
and their correlation functions that the anomaly kick, even with a tiny amplitude, acts as 
an effective trigger for the formation of the DCC domains. 

Our manuscript is organized as follows. In Section II we will review relevant aspects of 
the dynamics of the axial anomaly interaction in relativistic nuclear collisions. In Section III 
we discuss various issues of our numerical calculations, in particular, the initial conditions for 
the chiral fields and the method used to solve the field equations. In Section IV we present 
and discuss our numerical results: the effect of the anomaly kick on domain formation, the 
space-time evolution of the vector and axial currents, the evolution of the pion correlation 
functions, and the pion density distributions. We summarize our results and make some 
suggestions for future experiments in Section V. 

II. EFFECT OF CHIRAL ANOMALY IN DOMAIN FORMATION IN DCC; 

ANOMALY KICK 

We review the discussion in Ref. [25], in which two of us raised the question of possible 
effects of the chiral anomaly in triggering DCC formation in heavy ion collisions. There the 
influence of the strong electromagnetic fields of colliding ions on the isospin orientations of 
fields was discussed in the framework of the linear sigma model. This is a natural question 
to ask because electromagnetism breaks isospin. However, an explicit one-loop computation 
of the effective potential in the linear sigma model reveals that the isospin orientation of the 
ground state is not affected by the uniform background electromagnetic fields. Moreover, 
it can be argued on the basis of symmetry considerations that this result remains valid for 
any spatially and temporally varying background electromagnetic fields. Therefore, it was 
apparent that no effect is expected at the one-loop level. 

However, the above argument contains an interesting loophole. It is the presence of the 
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Wess-Zumino term which represents the effect of the chiral anomaly in the effective theory 
of QCD. It induces the it — > 27 coupling in the Lagrangian of the linear sigma model as 

£ = \d^i- \{\$\ 2 - fl) 2 + Ha + ^E-Htts, (1) 

where <fi = (er, 7F). The ground state of the chiral field therefore becomes sensitive to the 
background electromagnetic fields through the anomaly term. 

Then, what is the effect of the anomaly term on the chiral orientation of the pion fields, 
and in particular on the formation of DCC domains? To answer the question we first note 
that the equation of motion of the 1T3 field is altered to 

\2 ^ a 



II3 = V 2 7T 3 - A - f>) 7T 3 + —E-H, (2) 

where II3 denotes the conjugate field of 7r 3 . Then, II3 is affected by the strong electromagnetic 
fields of heavy ions during a collision by the amount 

a /■-»-. 

AU 3 = — E-Hdt, (3) 



where we have used the fact that the displacement of is negligible during the collision as its 
duration is quite short at RHIC and LHC P5| |. If we use the point-nucleus approximation, 



E-H can be easily evaluated by computing the Lienard-Wiechert potential. One obtains 

> - 2Z 2 e 2 ^ 2 

e - b = JWm1M (f - l) - (4) 



2 



\ 



#Tt*) ! + (?±t|) , (5) 



where we have assumed a collision of identical nuclei with mass M and charge Z, and 
7 denotes the Lorentz factor in the center of mass frame, and L = bx Mv with impact 
parameter b. 

The physical picture resulting from the above computation is as follows. The colliding 
two heavy ions have Lorentz contracted electromagnetic fields which produce a non-vanishing 
E-H for scatterings with angular momenta higher than s-wave. During the collisions the 
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interaction region of space-time is affected by the anomaly term prior to the DCC formation. 
Since the collision time scale R/j is much shorter than that of the DCC domains that may be 
formed in relativistic collisions, it is conceivable that it primarily affects the initial condition 
of the chiral fields. To visualize this effect of the anomaly term, we have introduced a term 
"anomaly kick" and described the effects imparted by the electromagnetic fields on the initial 
field configuration. 

The anomaly kick has a distinct characteristic feature which should be important in 
detecting its experimental signature. The II3 field is kicked toward opposite directions in 
isospin space in the upper and lower half-spaces off the scattering plane, as seen in (El). 



A rough estimate in Ref. | 25| , however, indicates that the magnitude of the kick is small, 
Z\n 3 ~ O.lrn" 2 . On the other hand, it is coherent over nuclear dimensions in each half-space. 



Therefore, while it was argued in Ref. [£5|] that the anomaly kick might produce interesting 
effects, it was difficult to draw definite conclusions on whether the anomaly effect enhances 
the domain formation of DCC. 

In this paper we implement the kick to the conjugate field of the neutral pion field as 
an initial condition of the numerical simulation code of DCC domain formation using the 
linear sigma model [IU|. 



We now elaborate the above estimate of AU3 by taking the finite size of the nucleus into 
account. While our subsequent treatment will not rely on any details of the spatial and 
temporal dependence of ^n 3 , it is critical for a reliable estimate of the magnitude of the 
kick induced by the chiral anomaly. Also, thus estimate will play a role in our discussions 
on the signature of the anomaly effect in DCC, which will be given in Sec. V. 

For the purpose of taking account of the finiteness of the nuclei, we consider a sphere 
with radius Ra inside which total charge Z is uniformly distributed, as an abstract of each 
colliding nucleus. We assume that two of such spheres are colliding at Lorentz factor 7 and 
impact parameter b in the center of mass frame, where each nucleus appears to be a sphere 
Lorentz contracted with the factor 7, or a pancake. We name the two pancake, pancake 1 
and pancake 2. In this case, E ■ B is given as follows: 



1) if r is inside pancake 1, i.e., R\ < R\, and inside the pancake of nucleus 2, i.e., 



- - 2Z 2 e 2 7 2 , A 

E ' B= wwk {fmL) ' (6) 



2) if r is inside pancake 1 and outside pancake 2, 



3) if r is outside pancake 1 and inside pancake 2, 



- - 2Z 2 e 2 -y 2 



4) if r is outside pancake 1 and outside pancake 2, 



- - 2Z 2 e 2 -y 2 

E - B =wmMm {f - L) - (9) 



Using Eqs. (^|) - @, we obtain ATl^ ~ 0.03m~ 2 in the central region of a non-central 
collision at RHIC, i.e., z ~ 0, h ~ r± ~ i?^ = i?(Au), r _L 6, and 7 ~ 100. The kick is 
substantially larger at LHC. 

III. NUMERICAL CALCULATION PRELIMINARIES 

A. Initial Condition and Constraints on Fields 

We have carried out numerical calculations on a lattice of 128 x 128, on which we have 
imposed periodic boundary conditions. The lattice is assumed to extend in the transverse 
directions perpendicular to the collision axis. In the longitudinal direction, the chiral fields 



are assumed to be boost invariant as in ||10|| . The lattice spacing is a. To include the initial 
correlation, we have adopted an initial correlation length, £ CO rr, which is generally larger 
than a. The initial fields and initial conjugate fields are set uniform within £ corr x £ corr cells 
and there is no inter-cell correlation. We have also carried out simulations with the quench 
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initial condition obtained from physically thermalized field configurations ]26[ and found 
that the result is insensitive to the details of the initial condition. 

The initial condition for the conjugate fields has been paid little attention to up to now. 
The initial conditions for both the chiral fields and the conjugate chiral fields, however, 
are expected to greatly affect the time evolution of the field configuration. Thus, we first 
consider the problem of the initial condition for the conjugate fields. 

We start with the following assumption on the chiral fields: the chiral fields can be 
treated statistically even in the case of the quench initial condition PJTof. This means that 
even if the initial field fluctuation is small, the fields have experienced some ergodic processes 
during the preceding phases of the heavy ion collision. 

The Hamiltonian density, 7i(x,t), of the linear sigma model above T c , is given by 

H(x,t) = t)\ 2 + i|vfe t)\ 2 + ^(|fet)| 2 + w 2 ) 2 - Ha, (10) 

We have replaced — v 2 by w 2 in order to take account of the chiral restoration in Eq. (TO). 
In the mean field approximation, w 2 and v 2 are related by flTC| , Pfl , 

w 2 = 3(,50j) + (50i)-t; 2 , (11) 

where Sfau (i = — 3) is the component of the fluctuation parallel to fa and 5<pi± is the 
orthogonal component. In fllQ|) , II is the conjugate field of the <fi field and should not be 
confused with the pion field. Since we are assuming that initially the system is in the 
symmetric phase and that the effective potential is convex, the motion of the chiral fields is 
obviously bounded. Accordingly, the virial theorem holds, 

^>^=(^Wj)) t > (12) 

where we have indicated that the average is a time average. Also note that there is no 
implied summation over multiple indices i in this subsection. Equation (|T2| ) leads to 

(Ui(x,t) 2 )t = (fa(x,t)(-V 2 + ^(x,t) 2 )fa(x,t)) t , (13) 

with 



fi(x,t) 2 = X(\$(x,t)\ 2 + w 2 ). (14) 

We have neglected a term proportional to H in (|I~3|). 

In numerical calculation, the Laplacian is discretized, and this affects the virial theorem, 
(|l~3|). In both of the methods used in our numerical calculations we shall discuss in the next 
subsection, the second derivative in the spatial directions is discretized by 

^72,/^ f(x + an i ) + f(x-an i )-2f(x) 

VJ(x) = , (15) 

Lb 

where is the unit vector in the i-th direction. For the correlation length £ co „ = 2a, we 
find 



i (x,t)V 2 i (x,t)) 4 = -^(0i(f,t) 2 ) t , (16) 

(X 



where D is the spatial dimension, and the virial theorem becomes 

( D 



^(<Mx,t) 2 ) t . (17) 

We assume that the field fluctuations at the initial time also satisfy the spatial analog of 
(0), 

(n 4 (f,t) 2 ) ccll W ^(<Mf,t) 2 ) C ell, (18) 
cr 

where the average is taken over uncorrelated cells and events. We use the value £ COTT = 2a 
in our numerical simulations throughout this paper. 



B. Calculational Procedure 

We have solved the equations of motion with the initial condition discussed in the pre- 
vious subsection. In numerical calculation, the equations of motion are discretized. It is 
well-known that naively discretized partial differential equations can be numerically highly 
unstable, even if the underlying original continuum version of the equations is stable. On 
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the other hand, what is essential in the time evolution of DCC domains is the amplifica- 
tion of low momentum modes, or in other words, the physical instability. Thus, we have 
to satisfy the following two apparently contradictory requirements in solving the equations 
of motion on discretized lattice at the same time: eliminate unphysical instabilities, but 
retain the physical ones. In order to achieve this goal, algorithms such as the Lax method 
and leap-frog method are often adopted. However, such algorithms introduce quite often 
numerical viscosity. As a result, unwanted spurious suppression of fluctuations of short wave 
lengths occurs. This does not cause a practical problem provided that the lattice spacing 
is small enough compared to the typical scale of the spatial variation of the fields. In the 
simulation of the formation of DCC domains, however, the behavior matters because the 
fluctuations of short wave lengths are responsible for the change in the effective potential flO 



and they affect the phase transition. Therefore, it is crucially important to keep short-range 
fluctuations to discuss the time evolution of DCC domains. 

In order to satisfy these requirements, we have used the following two methods, (i) 
the first order Adams-Bashforth method [^7j and (ii) the staggered leapfrog method for 
the second order term and the Crank-Nicholson method for the first order term in time 
28fl . Since (ii) is described in [F2^], we briefly explain (i) here. Detailed discussions on 



the numerical calculations will be presented elsewhere pq| . In both methods, if the lattice 



spacing is the same for all the spatial directions, the Laplacian is discretized as 

A/(zj, yj, T n ) -> (f(x i+1 , yj, r n ) + yj, r n ) 

+f(xi,y j+1 ,T n ) + f(xi,yj-i,T n ) - Af(x l ,y j ,T n ))/a 2 , (19) 



where r is the proper time defined by y/t 2 — z 2 . Note that only the two dimensional Lapla- 
cian appears in the equations of motion since we have assumed the longitudinal boost in- 



variance 



T0|.|T3[. The difference is in how to carry out the time integration. In the Adams- 
Bashforth method, the increment of /, f(xi, yj, r n+ i) — f(xi, yj, r ra ), is given by |2lJ 

1 5 

f(x h yj, T n+1 ) - f(x h yj, T n ) = Ar(f(x h yj, r n ) + - \jf{x h yj, r n ) + — XJ 2 f{xi, yj, r n ) 
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+ g V 3 J'(xi, yj, r n ) + — \7 4 f{xi, V j, T n ) + ■■■), (20) 

where At is the increment in r, f(xi, yj, r n ) is the proper time derivative of f(xi, yj, r n ), and 
\/ k is the backward difference operator defined by 

V 9(?i, yj, T n ) = g(xi, y d , r n ) - g(x h yj, r n -i), 
y k+l g(xi, yj, T n ) = v(V k g(%i, Vj, r n ) - V k g(x i} yj, r n _i)). (21) 

The time integrals of a(xi,yj,T n ), &(xi,yj,T n ), Tr(xi,yj,T n ), and n(xi,yj,T n ) are calculated 
at each lattice site with the help of (|20f) . In order to achieve satisfactory suppression of the 
unphysical instability, it is necessary to use sufficiently small At. For the typical parameters 
for the lattice spacing, the initial field fluctuation, and so on, which we shall specify shortly, 
we have found that At needs to be at least 1/1800 fm or smaller. We have found that 
the first order Adams-Bashforth method, which uses the first two terms in the series on the 
RHS of Eq. (|2Ti|), is good enough to suppress the unphysical instability. In the following 
calculations, we adopt the above value for At. We have also found that the required fineness 
of At in the method (ii) is similar and have verified that the two methods give similar results. 



IV. NUMERICAL CALCULATION RESULTS AND DISCUSSIONS 

A. Parameters 

In the following calculations, we shall use A = 19.97, v = 87.4 MeV, and H = (119 MeV) 3 . 
These values correspond to the pion mass m n = 135 MeV, the sigma mass m a = 600 MeV, 
and the pion decay constant /„- = 92.5 MeV at T = 0. We have assumed the longitudinal 
boost invariance |13[ at the initial time r , which is fixed at 1 fm, and we have fixed the 



lattice spacing a = 0.25 fm throughout the calculations. The initial <p % and <p l |3(| fields 
in each correlated cell are randomly distributed according to a Gaussian form. In order to 
take account of the finiteness of the initial hot system in the transverse directions, we have 
adopted the following Gaussian noise parameters: 

12 



(a) = (1 - f(r))f w , 
fa) = 0, 
(* 2 ) - = fa 2 > = S 2 f(r), 
(a) = fa) = fa) = 0, 
fa) = sgn(y) a n mlf(r), 

(- 2 ) = fa) = ^/(r), (22) 

where r is the distance from the origin to the center of a correlated cell and <5q is a constant. 
In the following calculations, we use 5q = v 2 /16, which corresponds to the quench scenario 



defined in [ 10 1 . In relating (a 2 ) and fa 2 ) to (a 2 ) and (n 2 ), we have taken advantage of the 



virial theorem (fTp) . f(r) is an interpolation function defined by 



/fa 



T Rq 

exp — - — + 1 



(23) 



and describes the boundary condition. R is the radius of the initially excited region where 
fluctuations of the classical chiral fields exist and the mean fields are different from their 
values in the vacuum. Outside this region, the chiral fields take the vacuum configuration, 
(/•n-,0). T is the thickness of the transient region. The results presented in this paper 
are obtained with T = 0.5 fm. a n has been introduced to take account of the effect of 
the anomaly kick discussed in Section [nj. a n is, in principle, dependent on the transverse 
coordinate. We, however, assume that it is a constant in order to concentrate on the effect 
of the kick to the time evolution of the system. Thus, the initial 7r 3 field consists of two 
parts, a randomly fluctuating part and a part with global coherence. In Eq. fl22|) , we have 
defined y direction to be perpendicular to the scattering plane. Following the estimate in 
Sec. II, we shall take a n = 0.1 in the calculation, which is of similar order but a bit larger 
than the one expected at RHIC, but much smaller than the one expected at LHC We note 
that this value is quite a moderate one as we see in the next subsection. 
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B. Results I — 



Effect of the Kick 



In Figs. 1 and 2, we show an example of initial tt 3 and 7r 2 field configurations, respectively. 
The x axis is parallel to the scattering plane, and perpendicular to the collision axis which is 
taken along the z direction. In this and the following calculations, we assume that the system 
extends infinitely in the transverse directions, i.e., R = oo, unless otherwise specified. 
With a n = 0.1 it is almost impossible to recognize the embedded coherent kick, (irs) = 
sgn(y)a n m^ in Eq. (|2"2]), in the configuration of the 773 field. In Fig. 3, we show another 
initial configuration of 7T 3 , where the value of the kick has been artificially increased by a 
factor of five to a n = 0.5. It is only for illustrative purpose; if the value of a n were that 
large, it would be possible to distinguish the initial configurations of 7r 2 and 773. However, 
such a large value of a n is not what is suggested by the real world at least at RHIC, and the 
value we use in our simulations is a n = 0.1. 

At a glance of Figs. 1 and 2, it appears that the effect of the anomaly kick may be too tiny 



to be observable. However, two of us argued in Ref. |25| that it might produce non-negligible 
effects on the evolution of the system because it is coherent over the nuclear dimensions. 
Moreover, the equations of motion for the chiral fields are coupled and highly non-linear, 
and therefore developing the physical intuition out of them would be quite difficult. Thus, 
we have decided to simulate the system numerically to uncover the time evolution of the 7r 2 
and 7T3 fields with the initial conjugate field configurations with the effect of the anomaly 
kick. It appears to the authors that it is the only tractable way of answering to the question 
of how the anomaly kick is effective in triggering DCC domains. 

The result we obtained is quite an unexpected one; the effect turns out to be significant. 
We show the proper time evolution of 113 and 7r 2 , respectively, in Figs. 4 and 5 every 3 fm 
from To in an event. Fig. 4 indicates that the small kick results in very prominent asymmetry 
of the 773 field between the upper and lower half-spaces. On the other hand, as is shown in 
Fig. 5, the 7r 2 field does not show any upper- lower space asymmetry in the course of the 
time evolution. We have also carried out calculation with a„ = and have found that the 
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time evolution of the 112 field is not affected by changing the value of a n Thus, the 

motion of different isospin components of pion fields effectively decouples from each other, 
and the anomaly kick affects only the time evolution of 1T3. We will explain the reason for 
this effective decoupling below. 

In addition to the existence of the upper-lower asymmetry, we find that the coherent 
structure which produces this asymmetry oscillates in time. The medium scale structure 
which both figures have in common corresponds to what is usually referred to as the DCC. 
Figs. 4 and 5 tell us that the anomaly kick induces a spatial structure in the field 
configuration which is larger than the ordinary DCC domains. This can also be observed in 
Figs. 6 and 7, in which we have plotted the Fourier power of 7r 3 and 7r 2 at r = r +3 fm = 4 fm, 
respectively. The scale is arbitrary, but the same scale has been used for Figs. 6 and 7. Note 
that this value of r corresponds to the second proper time in Figs. 4 and 5. Both figures 
show a strong enhancement of low momentum modes, but in Fig. 6 we observe further 
sharp enhancement at very low momenta. The former is due to the amplification of low 
momentum modes in the formation of normal DCC domains first discussed in Ref. ||, and 
the latter is due to the asymmetric coherent collective oscillation of the 7r 3 field induced by 
the anomaly kick. 

How can we understand such large effect in the time evolution of the 113 field caused 
by such a small kick? The reason is twofold. First, it is because the different field modes 
are effectively decoupled in the quench case. The equations of motion for the fields are 
non-linear, but the interaction term in the action is proportional to (0 l ) 2 (0- J ) 2 . Thus, it is 
expected that in the quench case, where the fluctuation of the fields is small, the approach 
toward the equipartition is slow. This is a simple analog of the fact that the relaxation time 
is longer at lower temperatures in the kinetic theory. Accordingly, low momentum modes 
behave as if they were almost independent oscillators undisturbed by the rest of the system. 
Also, in the quench scenario, the fields are initially concentrated at the local maximum of 
the effective potential with little fluctuation, and so even a small kick is expected to be 
effective in determining the subsequent motion of the fields if it is coherent in space. 
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The second reason is the weak coupling of Nambu-Goldstone modes at low momenta. 
When there exists expansion, the field configurations approach to the chiral circle quickly 
and gets concentrated around the chiral circle after a very short time PJT0|j. As is well- 
known, if the chiral fields are constrained on the chiral circle, the effective coupling among 
the pion fields becomes of the form of ~ (Tr-d^n) 2 . As a result, modes with small momenta 
are isolated from the rest of the system, and the effect of the coherent kick is not much 
affected by the presence of fluctuations of large momenta. 

C. Results II — Currents 

In the previous subsection, we have observed that the anomaly kick greatly affects the 
time evolution of the chiral fields. However, the field strengths themselves are not physical 
observables. The physical observables are currents, particle numbers, particle distributions, 
and so forth. In this subsection, we consider the behavior of physical currents. 

The 0(4) sigma model possesses an SU(2) x SU(2) symmetry in the limit of H —>■ 0. 
Correspondingly, two current densities, the vector current density V*(x, t) and axial current 
density A l (x, t), where i and ji are isospin and Lorentz indices, respectively, are defined as 

OH: 

V;(x,t)=e l ^(x,t)d tl 7c k (x,t), 

A^(x, t) = 7r* (s, t)d^a(x, t) - a(x, t)d^\x, t). (24) 
They satisfy the CVC and PCAC, respectively, 

d»v;(x,t) = o, 

d"Ai(x,t) = Hir\x,t). (25) 

In the following, we calculate averaged charge densities, i.e., the average of the zeroth com- 
ponent of the current densities. The average is taken separately over the upper and lower 
half-spaces with respect to the scattering plane. Average over events is taken also to reduce 
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fluctuations. Since the anomaly kick affects only the 713 field, its effect is expected to ap- 
pear only in the averages of Vq, V 2 , and Aq. Other charges are not expected to show any 
upper-lower asymmetries. 

In Fig. 8, we show (^l^upper, (^4o)u PP er, and (Al)\ ower , as a function of the proper time. 
The average is taken over 10 events in Figs. 8 - 10. We observe that Aq shows distinct 
upper- lower asymmetry. On the other hand, Aq does not show it. This, however, does not 
mean that Aq does not show domain structure. As is shown in Ref . , the low momentum 
modes of all the axial charges get amplified without the anomaly kick in the course of 
the time evolution from the quench initial condition. The reason for the lack of upper- 
lower asymmetry in the variable Aq is the absence of such an asymmetry in the initial field 
configuration. In Fig. 9, (V r 1 ) U p per and (lo 3 ) upper are shown. The vector charges do not show 
any asymmetry. Since the sign of the kick is fixed in our calculation, the expectation values 
of Aq in the upper and lower half-spaces do not vanish after taking average over events. 

The behaviors of (Ag) uppcr , (^4jj) upper ; and (Al)\ owcr are better understood by comparing 
them with the behaviors of (7ri) upp er, (Supper, and (7r 3 )i owcr , which are shown in Fig. 10. 
By comparing Fig. 8 and Fig. 10, we find that the extrema of (Al) npper and (Al)\ ower 
approximately correspond to the zero points of (7r 3 ) upper and (vr 3 )i owcr . This is what is 
expected from PCAC, which says that 

d°(Al(x, £)) = H(ir 3 (x, t)) + boundary terms. (26) 

An alternative and consistent interpretation of the above feature is that the time evolutions 
of (7r 3 ) upper and (7r 3 )i OW e r are approximated by harmonic oscillations as seen from Fig. 10, and 
thus at the zero points (Supper and (7r 3 )iower take extreme values. On the other hand, the 
sigma field does not receive any kick and, accordingly, does not show upper-lower asymmetry. 
Also, the oscillation of the sigma field is quickly damped and the sigma field approaches /„. 
T0| , |34|| . From this observation and the definition of the axial currents, Eq. it is thus 



concluded that the collective oscillation of the ir 3 field induced by the kick is responsible for 
the upper-lower asymmetry and the oscillatory behavior of (Aq). 
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The non-vanishing, upper-lower asymmetric, and oscillatory feature of (A 3 ) survives even 
if the transverse extent of the system is finite. In Fig. 11, we show (Al) upper , (Al) upper , and 
(^o) lower, for the case with Rq = 5 fm. The average is taken over all grid points with r < 5 
fm and over 100 events. In this case, the amplitudes of the oscillations of (v4o) U pper and 
(^4o)iowor are smaller than in the Ro = oo case. This is due to the reduction of the strengths 
of the pion fields caused by the transverse expansion. 



D. Results III — Correlation Functions 



In section [IV B| , we studied the time evolution of the chiral fields and their coherence in 



coordinate space visually by plotting the field strengths of the chiral fields. The best way of 
quantifying the coherence of fields is to calculate correlation functions. In this subsection, 
we study the correlation functions of the chiral fields and charges. 
First, we define the correlation function Aij(f, t) by 

Aij(r, t) = — J TTi(x, t)iXj(x + r, t)d z x — — J iti(x, t)d 3 x J 7tj{x\ t)d 3 x', (27) 

where V is the volume of the system. Since the fields 7Tj form the basis of a triplet represen- 
tation of SU(2), Aij(f,t) can be decomposed as |35 



Aij (f, t) = SijS(f, t) + e m V k (f, t) + T^r, t), (28) 

where S(f,t), Vk(r,t), and %j(f,t) are given by 

S{r,t) = ^A kk {f,t), 
Vk{r,t) = -£ ijk Aij{r,t), 

T^r, t) = ^{A i3 {r, t) + A^r, £)) - <%<S(f, t). (29) 

If Aij(f, t) = Aji(—r,t) holds, Vfc(f, t) vanishes. 

Next, let us assume a solution of the equation of motion for the 7Tj field, 7Tj(x, t), with 
given initial condition at t = to- Suppose the sign of the initial conditions for the (pi and fa 



fields is reversed, 4>i(x,t ) — > —(f>i(x,t ), 0j(5f, t ) — > —(pi(x,t ). Then, the solution for the 
7Tj field in this case is just —7ti(x,t). This is because the interaction term in the equation of 
motion for 7Tj is proportional to (a 2 + irf + 7r| + 7rf)7Tj. Thus, if the probability distribution 
P(TTk{x),7Tk{x)) of the initial values of the fields and -kk satisfies 

P(lT k (x), 7Tfe(x)) = P(-7T k (x), -7T fe (f)), (30) 

the event averages of all the off-diagonal correlation functions (Ai 3 -(f, £)) vanish, and hence 
(%j(f,t)) = for % ^ j. In the case we are dealing with, only the 713 field is kicked, and for 
the 7Ti and 7r 2 fields, the relationship (|30|) still holds. Accordingly, (7^(r, t)) = for i 7^ j. 
We have indeed confirmed this in our numerical calculation. Hence, in the following, we 
shall only consider the diagonal correlation functions. 

First, in Fig. 12, we show An(r,r) at three different proper times. We have used the 
default parameters and the result has been obtained by averaging over 5 events. The average 
over the azimuthal angle of r has been performed. We observe that the correlation function 
changes rather substantially as a function of the proper time. The pion fields oscillate and 
both the coherent part and short distance fluctuation change in time. Due to the energy 
conservation, when the amplitude of the coherent motion is at the maximum, the local 
fluctuation is at the minimum, r = 5 fm approximately corresponds to this time and in 
a two dimensional plot domain structure is clearly observed without fluctuations of much 
shorter scale. 

One interesting feature of An is the oscillatory behavior. This was already observed in 



Ref. [|K|. It has been sometimes interpreted as a sign of the shrinkage of the domains or 
occurrence of anti-correlation. However, this is not the correct interpretation. The behavior 
is actually due to the pseudo-periodicity in the distribution of the pion field strengths. This 
can be seen, for instance, in Fig. 5. The peaks of the pion field strength are not distributed 
randomly. The reason for this is traced back to the mechanism responsible for the formation 
of DCC domains, i.e., the amplification of low momentum pion modes. The wavelengths 
of the low momentum modes characterize not only the size of each domain but also the 
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distribution or separation of domains. In the mean field theory, modes with momentum 
less than a certain cutoff value are amplified at each time and not only specific modes are 
amplified. This is the reason why the distribution of the peaks is not utterly regular. We 
note that a similar oscillatory behavior in correlation functions of systems without regular 
lattice structure (such as liquids) is also known in condensed matter physics P5| . DCC 
formation is often compared to the spinodal decomposition observed in glass and metal 



37|,|38| . It is known, however, that in the case of the spinodal decomposition only modes 
in a small range of momenta are selectively amplified with large amplification factors. As a 
result, the distribution of domains created by the spinodal decomposition generally shows 
more periodic structure. 

Next, we show A 33 (r, r) at r = 7 and 9 fm in Fig. 13. As we have seen that the 
time evolution of the 7r 3 field is asymmetric between the upper and lower half-spaces, we 
distinguish the same side and different side correlation functions, A 33 ~(r,r) and A 33 ~(r, r). 
They are defined as follows: 

^33 + ( r > r ) = - (7r 3 (r)))(7r 3 (j,r) - (7r 3 (r))))y sgny(i) sgny(j) > 0, 

4b>,t) = ((7r 3 (i,r) - (7r 3 (r)))(7r 3 (j, r) - (7r 3 (r)»> - sgn y(i) sgny(j) < 0, (31) 

where % and j are grid points, (• • is the average over pairs of those grid points, % and j, 
such as the distance between them is r and the product of the signs of their ^/-coordinates 
is either positive (++) or negative (H — ), and (vr 3 (r)) is the average of the strength of the 
7r 3 field at r over the whole space. We show the average of Eq. ( [31]) over 5 events in Fig. 
13. From Fig. 10, we see the amplitude of the collective motion of the ir 3 field is about 
the maximum at r = 7 fm, and, accordingly, the upper-lower asymmetry is largest at that 
time. Correspondingly, we can clearly see the asymmetry in the correlation function, A 33 at 
t = 7 fm. On the contrary, at r = 9 fm, there is little asymmetry observed in v4 33 . This is 
because, as is seen from Fig. 10, the amplitude of the coherent oscillation almost vanishes 
at this time. For comparison, we plot the same side and different side correlations of the 7Ti 
field in Fig. 14 at r = 7 fm. As in the previous figures, an average over 5 events has been 
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taken. There is little difference between the two correlation functions. We have confirmed 
that this small amount of difference is solely due to statistical fluctuations by carrying out 
other sets of statistically independent simulations. 

We have further calculated the correlation functions of the axial and vector charges. 
We have defined the same side and different side correlation functions of the axial charges, 
Bki(r,r), and the vector charges, C&j(r, r), by 

B kl (r,r) = ((A k Q (i,r) - (A k (r)))(A l (j,r) - (A l (r)))) kl , 

C kl (r,r) = <(Kf(z,r) - (V k (r)))(V^r) - (^(r)»>„, (32) 

where the notations are the the same as before and the average is taken over 5 events. In 
Fig. 15, we show the same side and different side correlation functions of Aq at r = 9 fm. 
We note that this proper time corresponds to when the asymmetry in A^ is approximately 
at maximum as seen in Fig. 8. In Fig. 15, the asymmetry between the same side and 
different side correlation functions is observed. However, it should be noted that the degree 
of coherence is less than that of the field-field correlation functions. In Fig. 16, same side 
correlation functions of Aq, Vq, and V 3 are plotted at r = 9 fm. This indicates that A\ 
shows some coherence, but that neither Vq nor Vq 3 shows coherent behavior [33 . 



E. Results IV — Pion Density Distributions 

It has often been argued that the pion density in coordinate space is proportional to the 
square of the pion field strength. This is actually the basis for the expression for the proba- 
bility to have a DCC domain in which the fraction of neutral pions, n n o / (n^o + n n + + n n - ) , 
takes the value R |5|-[7]|. 

P(R) = -4=. (33) 

In this formula, the pion fields are assumed to be static, tt^ = 0. In reality, however, 
the pion fields are not static and also include substantial fluctuations. When the pion 
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fields are not static, the pion densities depend also on the conjugate fields. This is just an 
analog of the one dimensional harmonic oscillator case; both the coordinate variable and 
its canonical conjugate (momentum) variable contribute to the energy, i.e., the number of 
quanta. Moreover, when the pion fields are fluctuating even classically, the definition of a 
domain itself is not unambiguous. However, in the following we shall calculate the local 
pion densities in coordinate space as a function of proper time and discuss the typical size 
of clusters defined by the distribution of local pion densities thus obtained. 

Quantum field theory is often formulated in momentum space. Creation and annihilation 
operators are usually defined for momentum eigenstates. It is because the Hamiltonian 
for free fields reduces to the sum of those of harmonic oscillators in momentum space. 
However, as we have discussed above, what we need here is the pion density distribution 
in coordinate space. For that purpose, we first define the pion creation operators and 
annihilation operators in coordinate space as the Fourier transforms of those in momentum 
space. For example, the creation operators for n + , 7i~ , and ir° at x are given, respectively, 
by 

- (/P® - w"®) d3p ' 

- Wp/ e ~" (V?"-*® - 8 ;if ,(p_) ) d3p ' 

4.M = p^/— (V?><0 " A, (34) 

where u p = \Jm\ + p 2 , iti{p) and 7Tj(p) are, respectively, the Fourier transforms of 7Tj(x) and 
7Tj(x). The charge eigenstate fields n + (x), n~(x), and n°(x) are defined as 

7Ti(f) ±m 2 {x) 

71 {x) = T2 — ' 

7T°(£) =7T 3 (f). (35) 

We define the local density operators for the Hi field, n ni (x) as 

n v .{x) = al.(x)a v .(x). (36) 

We note that since 
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kHp) ± ^(P), (37) 

in general, 

n n +(x) 7^ n w - (x). (38) 

Thus, in principle, the electric charge density can fluctuate and can take non-zero values 
locally. The creation operators and density operators satisfy the following commutation 
relation: 

[n Ki (x), al.(x')} = 8 i:j 8(x- x')al.(x'). (39) 

This justifies our definition of the creation, annihilation, and density operators in coordinate 
space. These creation operators do not create momentum or energy eigenstates, since the 
position is specified. However, Parsival's relation is satisfied, 

n^(x)d 3 x = [n Wi {p)d 3 p, (40) 



and, accordingly, the integrations over coordinate space and momentum space give the same 
total number. 

In this subsection, we have not yet specified the physical state of the system. In order 
to extract the coordinate space distribution of each pion species and find the typical size of 
DCC domains in the original sense, we first rewrite follows: 



UirAX 



1 r 2 
(P(x - f>i(x') + iQ(x - x r )TXi{x')) d 3 x' 



(41) 



where P(x) and Q(x) are functions defined by 



P(x) = 7 ^UlJ- p e^d 3 p, 



1 

(2~) 

Q(x) = -Lj / -U^dV (42) 



The non-local feature of n ni is obvious in Eq. fl4~l"|). We now assume that the system is 
in a coherent state | c) and the operators n(x) and ii(x) in a(x) operating on | c) can be 
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replaced by those values obtained in the classical simulation, 7r c i(x) and tc c i(x). Note that 
this procedure needs to be modified if the system is not in a coherent state |2(J. This 
possibility will be considered elsewhere. 

In Figs. 17 and 18, we show the distributions of tt° and tt + density on the x-y plane, 
respectively, every 3 fm from r in an event. We have used the same initial field configuration 
as for Figs. 4 and 5. From Figs. 17 and 18, we can clearly see the formation of domains 
also in this definition. However, we recognize in these figures several features different from 
Figs. 4 and 5. First, the distributions of 7r° and 7r + densities fluctuate more than those 
of the 7T° and tt + fields themselves. This is because the definitions of tt° and 7r + density 
include, respectively, 7r° and 7r + fields, which are not so much amplified as tt° and tt + fields 
and noisier than those. Second, the distribution of 7r density does not show the upper-lower 
asymmetry. This is because the definition of n ni (x) is quadratic in tt^x) and 7Tj(x). Third, 
after the initial formation of the large scale structure, the pattern of the distribution remains 
almost unchanged except for the occurrence of slight diffusion. This tells us that domains 
in this definition are almost frozen at fixed positions on the x-y plane. The last point is 
also clearly seen in Fig. 19, where we have plotted the time evolution of the distribution of 
R(x) = n n o(x)/(n n +(x) +n 7T -(x) +n 7r o(x)) in the same event. 

In Fig. 20, we have plotted R, which is defined as 

R= J2 I ' n^{k)d?k I J(n w+ (k) +n n -(k) +n w o(k))d 3 k, (43) 

event event 

where n ni (k) is the pion density at momentum k, in a solid line. Ten events were used in 
the average. The dashed line is the same except that the momentum integration is limited 
to k < 250 MeV. Fig. 20 shows striking enhancement of neutral pions, especially in the 
low momentum region. Note that R = 0.375 corresponds to N n o : N 7T + : N n - = 1.2 : 1 : 1 
provided N n + = N w -. Also note that the distribution of R deviates from the celebrated 
form, P(R) = l/2yR, when 7ti(x) is not uniform, and moreover that even if 7T;(af) is spatially 
uniform, when 7fj(x) 7^ 0, R also depends on iii{x). Fig. 21 is the same as Fig. 20 except that 
no anomaly kick is exerted. In this case, R is, as is expected, approximately 1/3 regardless 
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of the cutoff in the momentum integral. 

We have found that the distribution of the excess 7T°'s is not isotropic in momentum 
space. They are distributed more in the y-direction than the x-direction. This reflects the 
Fourier power of the initial anomaly kick: in the Rq — > oo limit, the Fourier power of the 
kick, sgn(y)a n m^, is concentrated at p x = in the x-direction, while it has support at p y ^ 
in the ^/-direction. In addition, the effect of the anomaly kick is expected to be observed 
more in the central rapidity region than in the fragmentation region, because the kick is 
created when the two colliding nuclei are almost overlapping with each other. 

V. SUMMARY 

In summary, we have studied the physical influence of the axial anomaly on the formation 
mechanism of disoriented chiral condensate domains in relativistic heavy ion collisions. Our 
study is based on the framework of the linear sigma model with initial conditions correspond- 
ing to a "quench" modified by a small, spatially coherent displacement ( "anomaly kick" ) of 
the conjugate field of the neutral pion field caused by the action of the electromagnetic fields 
of the colliding nuclei. These initial conditions correspond to a somewhat schematic, but 
quantitatively realistic abstraction of conditions that could arise in semiperipheral collisions 
of two heavy nuclei at RHIC or LHC. 

We have found that the anomaly kick has practically no influence on the evolution of the 
charged pion fields but causes noticeable modifications to the dynamics of the neutral pion 
field. It produces a spatial asymmetry in the formation of DCC domains above and below 
the scattering plane and leads to a general enhancement of the 7r° domain structure. The 
up-down asymmetry is clearly expressed in the integrated I3 = component of the axial 
charge on both sides of the scattering plane. Since the signs of the pion fields or the axial 
charges are not directly connected with hadronic observables, this up-down asymmetry is 
difficult to establish experimentally. 

However, the overall enhancement of domain formation in the ^3 = direction is found to 
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lead to a relative increase in the neutral pion yield, raising the ratio R = N^o/N^ significantly 
above the value 1/3 dictated by isospin symmetry. This signature, which survives event 
averaging, should be experimentally detectable. Because it has a different characteristic 
than the isospin breaking due to the effects caused by charges of incident nuclei, which alters 
the 7r + /7r~ ratio, it can be discriminated from this well known effect. Also, the mechanism 
of isospin breaking discussed here can be distinguished from the pion mass splitting, which 
leads to a similar but smaller enhancement of R, because it depends on the impact parameter 
and vanishes for precisely central collisions. 

Finally, we have analyzed spatial correlation functions of the pion field. The oscillations 
already observed by others appear to be caused by long-range correlations in the positions 
of different DCC domains. These correlations correspond to a fluid-like spatial structure 
of the domain locations and are clearly visible in equal-time snapshots of the chiral order 
parameter. 

If the domain formation of DCC is in fact triggered by the anomaly effect discussed in 
this paper, there is a clear experimental signature which testifies for its origin; the formation 
rate of DCC should show a strong energy and atomic number dependences characteristic 
to the mechanism. The parameter a n which measures the magnitude of the anomaly kick 
scales as E and Z 2 , respectively, as center of mass energy E and the atomic number Z are 
varied. It is tempting to speculate that the formation rate of DCC somehow should reflect 
these dependences. 

Experimental searches for DCC formation have, so far, concentrated on central collisions 
between two heavy nuclei, because it is here where one expects the highest energy densities to 
occur over the largest spatial volumes. These two conditions are important, because matter 
needs to be heated above T c over a large volume if one wants to observe the characteristic 
effects of a phase transition. 

However, our results indicate that semi-peripheral collisions, i.e., collisions with an im- 
pact parameter b ~ Ra, may be even more favorable for DCC formation. Two reasons 
contribute to this phenomenon: (1) The expansion time scale is generally shorter in periph- 
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eral collisions than in central ones, and a three-dimensional expansion may occur more easily 
in non-central collisions, leading to the quench initial condition. (2) Collisions of heavy nu- 
clei with b ^ produce the "anomaly kick" , which favors the formation of large uniform 
domains of the quark condensate by inducing a collective motion in the direction of the 
neutral pion field. Therefore, we suggest that the search for DCC formation in relativistic 
heavy ion collisions at RHIC should not be confined to central collisions only. 
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FIGURE CAPTIONS 



Fig. 1: Example of initial 7r 3 field configuration. 
Fig. 2: Example of initial 7T2 field configuration. 

Fig. 3: Example of initial 7r 3 field configuration. The value of the kick has been artificially 

increased by a factor of five to a n = 0.5. 

Fig. 4: Proper time evolution of the 7r 3 field in an event. 

Fig. 5: Proper time evolution of the 7T2 field in an event. 

Fig. 6: Fourier power of 7r 3 at r = 4 fm in an event. 

Fig. 7: Fourier power of 7r 2 at r = 4 fm in an event. 

Fig. 8: (Al) uppcr , (Al) uppcr , and (-Aq) lower as a function of proper time. The average is taken 
over 10 events. 

Fig. 9: (V 1 ) U p pcr , (V 3 ) U ppcr as a function of proper time. The average is taken over 10 events. 
Fig. 10: (7ri) upper , (7r 3 ) upper , and (7r 3 )i ower as a function of proper time. The average is taken 
over 10 events. 

Fig. 11: (Al) uppcr , (Al) upper , and (A|)i ower as a function of proper time for the case with 
R = 5 fm. The average is taken over 100 events. 

Fig. 12: A u (r, r) at r = 1, 5, and 9 fm. The average is taken over 5 events. 

Fig. 13: A^(r, r) and A^(r, r) at r = 7 and 9 fm. The average is taken over 5 events. 

Fig. 14: Af-f(r, r) and Af^~(r, r) at r = 7 fm. The average is taken over 5 events. 

Fig. 15: i?33 + (r, r) and B^f(r,r) at r = 9 fm. The average is taken over 5 events. 

Fig. 16: B^{r,r), C{i + (r, r), and C^ 3 + (r,r) at r = 9 fm. The average is taken over 5 

events. 

Fig. 17: Proper time evolution of the distribution of n° density in an event. 
Fig. 18: Proper time evolution of the distribution of rt + density in an event. 
Fig. 19: Proper time evolution of the distribution of R(x). 

Fig. 20: Solid line is the proper time evolution of R averaged over 10 events. The dashed 
line is the same except that the momentum integration is limited to k < 250 MeV. 
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Fig. 21: Same as Fig. 20 except that no anomaly kick is exerted. 
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